Code
knitr::opts_chunk$set(
message = FALSE
)Supplementary Report S2. Report with ASV without OTU post-clustering. Terroir and farming practices drive arbuscular mycorrhizal fungal communities in French vineyard.
Code
#if (!require("devtools", quietly = TRUE)) {
# install.packages("devtools")
#}
#devtools::install_github("adrientaudiere/MiscMetabar")Code
library(MiscMetabar)
library(targets)
library(metacoder)
library(DT)
library(knitr)
library(ggplot2)
library(patchwork)
library(vegan)
library(FactoMineR)
library(factoextra)
library(DESeq2)
library(indicspecies)
library(sf)
library(adespatial)
library(ade4)
library(adegraphics)
library(kableExtra)
library(formattable)
library(microbiomeMarker)
library(ComplexUpset)Code
d <- tar_read(d_blast)Metadata modification
Code
d@sam_data$rank <- gsub("Weeding", "Grassing", d@sam_data$rank)
d@sam_data$inter_rank <- gsub("Weeding", "Grassing", d@sam_data$inter_rank)
d@sam_data$soil_practice <- gsub("Weeding", "Grassing", d@sam_data$soil_practice)
d@sam_data$terroir <- gsub("Perrier", "Vergèze", d@sam_data$terroir)Spores number
Code
d@sam_data$nb_spores <- rowMeans(cbind(as.numeric(d@sam_data$spores_rep1), as.numeric(d@sam_data$spores_rep2),as.numeric(d@sam_data$spores_rep3)),na.rm=TRUE)Practices
Code
d@sam_data$organic <- ifelse(!d@sam_data$practice %in%
c("Conventional", "Conversion"),
"Organic", "NonOrganic"
)Paired samples (roots, spores from the same samples)
Code
d@sam_data$paired_name <-
apply(as.matrix(d@sam_data[, c(5:16)]), 1, paste,
collapse = "--"
)
d@sam_data$paired_name[d@sam_data$terroir == "Vergèze"] <-
gsub("R", "", d@sam_data$ref_mycea[d@sam_data$terroir == "Vergèze"])Code
tib_sam_data <- as_tibble(d@sam_data) %>%
mutate(across(starts_with("Myc_freq"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_colonization"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_root"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_richness"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_abundance"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_richness"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_abundance"), as.numeric)) %>%
mutate(Myc_freq = rowMeans(select(., starts_with("Myc_freq")), na.rm = TRUE)) %>%
mutate(Myc_intensity_colonization = rowMeans(select(
., starts_with("Myc_intensity_colonization")
), na.rm = TRUE)) %>%
mutate(Myc_intensity_root = rowMeans(select(., starts_with(
"Myc_intensity_root"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_richness = rowMeans(select(., starts_with(
"Arbuscul_richness"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_abundance = rowMeans(select(., starts_with(
"Arbuscul_abundance"
)), na.rm = TRUE)) %>%
mutate(Vesicle_richness = rowMeans(
select(., starts_with(
"Vesicle_richness"
)),
na.rm = TRUE
)) %>%
mutate(Vesicle_abundance = rowMeans(
select(., starts_with(
"Vesicle_abundance"
)),
na.rm = TRUE
)) %>% tibble::column_to_rownames(var = "X")
sample_data(d) <- tib_sam_dataExclude Mycelium
Code
d <- clean_pq(subset_samples(d, compartment != "Mycelium"))
d_vs <- clean_pq(subset_samples(tar_read(d_vs), compartment != "Mycelium"))Summary of final dataset
Code
summary_plot_pq(d) +
labs(title = "Summary of final dataset") +
theme(plot.title = element_text(hjust = 0.5, size = 20, color = "#1e2b4c"))
Summary of bioinformatics pipeline
Global summary
Code
kable(tar_read(track_sequences_samples_clusters))| nb_sequences | nb_clusters | nb_samples | |
|---|---|---|---|
| Paired sequences | 8283716 | 28637 | 194 |
| Paired sequences without chimera | 7777304 | 9791 | 194 |
| Paired sequences without chimera and longer than 200bp | 7736170 | 7522 | 194 |
| ASV denoising | 7533009 | 7522 | 186 |
| OTU after vsearch reclustering at 97% | 7533009 | 1081 | 186 |
| OTU after blast filter without reclustering | 3819912 | 3793 | 182 |
| OTU after blast filter with reclustering | 3822304 | 219 | 182 |
SUBSETING AND FILTERING DATASET
Code
d_muco <- clean_pq(subset_taxa(d, Class_PR2 == "Mucoromycota"))
d_blast <- clean_pq(subset_taxa(tar_read(d_blast), Class_PR2 == "Mucoromycota"))Code
for(i in c(19:46)){
d_muco@sam_data[,i] <- as.numeric(unlist(d_muco@sam_data[,i]))
}Code
d_roots <- clean_pq(subset_samples(d_muco, compartment == "Roots"))
d_spores <- clean_pq(subset_samples(d_muco, compartment == "Spores"))Code
d_rs_merged_samples <- clean_pq(merge_samples2(
d_muco,
"paired_name",
default_fun =
function(x){MiscMetabar::diff_fct_diff_class(x, na.rm = T)}
))Code
d_rs <- clean_pq(subset_samples_pq(
d_rs_merged_samples,
sample_sums(d_rs_merged_samples) > 2000
))TAXONOMICAL ANALYSES
Code
# Figure S8
multitax_bar_pq(d_vs, "Supergroup_PR2",
"Subdivision_PR2", "Class_PR2",
nb_seq = TRUE
) +
ggtitle("Number of sequences (log10) including non-Mucoromycota")
Number of sequences
Code
# Figure S6
dir.create("output/krona")Warning in dir.create("output/krona"): 'output/krona' existe déjàCode
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk",
name = "OTU_post_clustered_allEuk"
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM",
name = "OTU_filtOnMaarjam_AM"
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM",
name = "OTU_filtOnMaarjamAndPr2_AM"
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk",
"output/krona/OTU_filtOnMaarjam_AM",
"output/krona/OTU_filtOnMaarjamAndPr2_AM"
),
"output/krona/Euk_to_AM_filtering_nb_seq.html"
)Number of OTU
Code
# Figure S7
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk_Nbotu",
name = "OTU_post_clustered_allEuk",
nb_seq = F
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
name = "OTU_filtOnMaarjam_AM",
nb_seq = F
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu",
name = "OTU_filtOnMaarjamAndPr2_AM",
nb_seq = F
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk_Nbotu",
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu"
),
"output/krona/Euk_to_AM_filtering_Nbotu.html"
)NUMBER OF SPORES ANALYSES
Code
psm_otu <- psmelt(as_binary_otu_table(d_spores)) |>
filter(Abundance > 0) |>
group_by(Sample) |>
summarise(
"Abundance" = sum(Abundance),
"region" = region[1],
"nb_spores" = nb_spores[1]
)Code
psm_res <- psmelt_samples_pq(d_spores) %>%
mutate(nb_spores_log10 = log10(nb_spores))Spores count <-> Alpha diversity
Code
# Fig S2
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_spores, rngseed = 221013)) %>% mutate(nb_spores_log10 = log10(nb_spores))
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Abundance, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_0, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_1, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_2, type = "non-parametric")
Spores count <-> Practice
Code
# Figure S5
ggstatsplot::ggbetweenstats(psm_res, practice, nb_spores_log10, type = "non-parametric") 
Spores count <-> Terroir
Code
# Figure S4
ggstatsplot::ggbetweenstats(
mutate(psm_res,
terroir = reorder(psm_res$terroir, psm_res$nb_spores)
),
terroir,
nb_spores_log10,
type = "non-parametric",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist"
) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
Spores count <-> Mycorrhization rate
Code
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_freq, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_intensity_colonization, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Arbuscul_abundance, type = "non-parametric")
MYCORRHIZATION MEASURE
Mycorrhization rate <-> Alpha diversity
F = Myc_freq
Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))
# Figure S3a
(ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_2, type = "non-parametrique"))
M = Myc_intensity_colonization
Code
# Figure S3b
(ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_2, type = "non-parametrique"))
A = Arbuscul_abundance
Code
# Figure S3c
(ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_2, type = "non-parametrique"))
Mycorrhization rate <-> terroir
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)
Mycorrhization rate <-> practice
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)
ECOLOGICAL ANALYSES : ALPHA DIVERSITY
Compartment
Code
psm_res <- psmelt_samples_pq(d_muco)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_muco, rngseed = 221013))Code
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_0,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_1,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_2,
type = "non-parametrique"
)
Code
ggplot(psm_res_rarefied, aes(y=Hill_0, x=compartment, color=practice))+
geom_point() +
geom_line(aes(group = paired_name)) +
facet_wrap(~terroir)
Terroir
Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))Code
# Figure 4
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")
Code
# Figure S10
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_0))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_1))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_2))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")
Code
# Fig_7a
p_accu_balanced_mod_terroir <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"terroir",
nperm = 999,
step = 300,
quantile_prob=0.9,
sample.size = 10000,
progress_bar = FALSE)Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyéCode
p_accu_balanced_mod_terroir +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod_terroir$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod_terroir$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,25000)) + theme(legend.position = "none" )Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).
Practice
Code
# Figure 5
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 0 (richness)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 1 (~Shannon)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 2 (~Simpson)") + plot_layout(axes = "collect")
Code
# Figure S11
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
)
Code
# Fig_7b
p_accu_balanced_mod <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"practice",
nperm = 999,
step = 300,
quantile_prob=0.90,
progress_bar = FALSE)Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyéCode
p_accu_balanced_mod +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,18000)) + theme(legend.position = "none" )Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).
ZOOM ON LOW-DIVERSITY SAMPLES
Code
otu_hill <- vegan::renyi(d_rs@otu_table, scale = c(0, 1, 2), hill = TRUE)
low_div_samples <-
rownames(otu_hill[otu_hill$`0` < 10 |
otu_hill$`1` < 2 |
otu_hill$`2` < 2, ])
low_div_samples[1] "PE13"
[2] "2020--Bordelais--Langoiran--Langoiran--Conventional chimique--Conventional---0.356503--44.713204--Conventional--Herbicide--Herbicide--Herbicide"Code
sum(otu_hill$`0` < 10)[1] 2Code
sum(otu_hill$`1` < 2)[1] 0Code
sum(otu_hill$`2` < 2)[1] 0Code
length(low_div_samples)[1] 2Code
low_div_samples_table <- tibble(cbind(as_tibble(d_rs@sam_data[low_div_samples, c("terroir", "practice", "rank", "inter_rank", "compartment")]), apply(otu_hill[low_div_samples, ], 2, round, digits = 2), sample_sums(d_rs)[low_div_samples]))Warning in class(x) <- tibble_class: La class(x) est une chaîne multiple
("tbl_df", "tbl", ...) ; le résultat ne sera plus un objet S4Code
low_div_samples_table$compartment <- ifelse(is.na(low_div_samples_table$compartment), "Spores + Roots", "Spores")
colnames(low_div_samples_table) <- c("Terroir", "Global practice", "Rank practice", "Inter rank practice", "Compartment", "Hill 0 (richness)", "Hill 1 (~Shannon)", "Hill 2 (~Simpson)", "nb_seq")
low_div_samples_table <- arrange(low_div_samples_table, as.numeric(`Hill 0 (richness)`))
low_div_samples_table <-
rbind(
low_div_samples_table,
c(
"",
"",
"",
"MEAN",
"(low -diversity samples)",
round(mean(low_div_samples_table$`Hill 0 (richness)`), 2),
round(mean(low_div_samples_table$`Hill 1 (~Shannon)`), 2),
round(mean(low_div_samples_table$`Hill 2 (~Simpson)`), 2),
round(mean(low_div_samples_table$nb_seq), 2)
),
c(
"",
"",
"",
"MEAN",
"(all samples)",
round(mean(otu_hill$`0`), 2),
round(mean(otu_hill$`1`), 2),
round(mean(otu_hill$`2`), 2),
round(mean(sample_sums(d_rs)), 2)
),
c(
"",
"",
"",
"MAX",
"(all samples)",
round(max(otu_hill$`0`), 2),
round(max(otu_hill$`1`), 2),
round(max(otu_hill$`2`), 2),
round(max(sample_sums(d_rs)), 2)
)
)
# Table 3
kbl(low_div_samples_table) |>
kable_classic(full_width = F, html_font = "Cambria")| Terroir | Global practice | Rank practice | Inter rank practice | Compartment | Hill 0 (richness) | Hill 1 (~Shannon) | Hill 2 (~Simpson) | nb_seq |
|---|---|---|---|---|---|---|---|---|
| Langoiran | Conventional | Herbicide | Herbicide | Spores + Roots | 6 | 3.05 | 2.11 | 12643 |
| Vergèze | Organic | Scratching | Ploughing | Spores | 7 | 4.86 | 3.95 | 5938 |
| MEAN | (low -diversity samples) | 6.5 | 3.96 | 3.03 | 9290.5 | |||
| MEAN | (all samples) | 101.28 | 20.09 | 10.46 | 48265.19 | |||
| MAX | (all samples) | 432 | 64.2 | 24.15 | 101104 |
ECOLOGICAL ANALYSES : BETA DIVERSITY
Compartment
Code
# Figure S9
p <- phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), "matrix"),
method = "robust.aitchison"
)),
color = "region",
shape = "compartment"
) +
geom_point(size = 3) +
stat_ellipse(inherit.aes = F, aes(x = DCA1, y = DCA2, linetype = compartment))Warning in phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), : `Ordination species/OTU/taxa
coordinate indices did not match `physeq` index names. Setting corresponding
coordinates to NULL.Code
p + geom_line(data = p$data, aes(group = paired_name))
Practice and terroir
Code
p <-
phyloseq::plot_ordination(d_rs,
ordinate(d_rs,
method = "NMDS",
distance = "bray"
),
color = "terroir",
shape = "practice"
) + facet_wrap("region") + scale_shape_manual(values=c(25,23,24))Square root transformation
Wisconsin double standardization
Run 0 stress 0.2199938
Run 1 stress 0.2203704
... Procrustes: rmse 0.03925079 max resid 0.2885775
Run 2 stress 0.2265037
Run 3 stress 0.2213033
Run 4 stress 0.2203289
... Procrustes: rmse 0.03921461 max resid 0.2889235
Run 5 stress 0.2201826
... Procrustes: rmse 0.07042354 max resid 0.5261175
Run 6 stress 0.242704
Run 7 stress 0.2304432
Run 8 stress 0.220005
... Procrustes: rmse 0.07747605 max resid 0.510776
Run 9 stress 0.2200585
... Procrustes: rmse 0.01712057 max resid 0.1278995
Run 10 stress 0.2209168
Run 11 stress 0.2241108
Run 12 stress 0.2203529
... Procrustes: rmse 0.0391036 max resid 0.2892524
Run 13 stress 0.2192479
... New best solution
... Procrustes: rmse 0.07702747 max resid 0.514501
Run 14 stress 0.2196043
... Procrustes: rmse 0.03466159 max resid 0.2841357
Run 15 stress 0.2214039
Run 16 stress 0.2216717
Run 17 stress 0.2447455
Run 18 stress 0.2183676
... New best solution
... Procrustes: rmse 0.03750409 max resid 0.2833134
Run 19 stress 0.2217418
Run 20 stress 0.2201252
*** Best solution was not repeated -- monoMDS stopping criteria:
1: no. of iterations >= maxit
19: stress ratio > sratmaxCode
p_nmds <- p + geom_point(
data = select(p$data, -c(region)),
fill = "grey40",
color = "grey20",
size = 1,
alpha = 0.5
) +
geom_point(size = 3, aes(fill=terroir)) +
ggtitle("NMDS on bray distance") +
scale_fill_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0")) +
scale_color_manual(values=rep("grey20",14))
# Figure 9
p_nmds
Code
# Figure 6a
upset_pq(d_rs, "practice", taxa_fill = "Family",
set_sizes=(
upset_set_size()
+ geom_text(aes(label=after_stat(count)), hjust=-0.3, color="white", stat='count')
))
Code
ggvenn_pq(d_rs, "practice")
Code
# Figure 6b
upset_pq(d_rs, "terroir", taxa_fill = "Family", min_size = 2, height_ratio = 0.6
,set_sizes=(
upset_set_size()
+ geom_label(aes(label=after_stat(count)), hjust=-0.3, size=2.5, stat='count')
))Warning: Removed 233 rows containing non-finite outside the scale range
(`stat_count()`).
Global beta-div analysis
Spatial data
Code
library(geodist)
lat_lon <- as_tibble(d_rs@sam_data[, c("lat", "long")])
lat_lon$lat <- as.numeric(lat_lon$lat)
lat_lon$long <- as.numeric(lat_lon$long)
dist_spatial_meter <- as.dist(geodist(lat_lon, measure = "geodesic"),
upper = FALSE
)Code
lon_lat_rs <- d_rs@sam_data[, c("long", "lat")]
lon_lat_rs$long <- as.numeric(lon_lat_rs$long)
lon_lat_rs$lat <- as.numeric(lon_lat_rs$lat)
MEM <- dbmem(lon_lat_rs, MEM.autocor = "non-null")Code
test_MEM <- moran.randtest(MEM, nrepet = 999)
test_MEM$pvalue_adjust <- p.adjust(test_MEM$pvalue, method = "BH")Code
d_rs@sam_data$MEM_1 <- MEM[, 1]
d_rs@sam_data$MEM_2 <- MEM[, 2]
res_ado_spatial_robAit <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_robAit_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
res_ado_spatial_bray <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "bray"
)
res_ado_spatial_bray_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "bray"
)
df <- data.frame(res_ado_spatial_bray)
df$names <-
factor(
rownames(df),
levels = c(
"sample_size",
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p1 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray")
df <- data.frame(res_ado_spatial_robAit)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p2 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison")
df <- data.frame(res_ado_spatial_bray_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p3 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray after rarefaction")
df <- data.frame(res_ado_spatial_robAit_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p4 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison after rarefaction")
(p1 + ylim(0, 0.12) + p2 + ylim(0, 0.12)) /
(p3 + ylim(0, 0.12) + p4 + ylim(0, 0.12))Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.23494 0.018072 1.4539 0.1618
Residuals 61 0.75823 0.012430 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.01739 0.0086972 0.8898 0.4152
Residuals 72 0.70378 0.0097748 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.23494 0.018072 1.4539 0.1618
Residuals 61 0.75823 0.012430 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.01739 0.0086972 0.8898 0.4152
Residuals 72 0.70378 0.0097748 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 2471.66 190.127 13.793 1.288e-13 ***
Residuals 61 840.86 13.785
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 66.0 32.985 0.4173 0.6604
Residuals 72 5691.6 79.050 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 2471.66 190.127 13.793 1.288e-13 ***
Residuals 61 840.86 13.785
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 66.0 32.985 0.4173 0.6604
Residuals 72 5691.6 79.050 Soil properties
Code
soil_prop <- as_tibble(d_rs@sam_data) |>
select(paired_name, practice, organic, terroir, Coarse_sand:CEC)
soil_prop_num <- soil_prop |>
select(-all_of(c("practice", "organic", "terroir", "CaO", "Total_sand", "Total_filt"))) |>
tibble::column_to_rownames("paired_name") |>
mutate(across(everything(), as.numeric)) |>
tidyr::drop_na()Code
pca_test_res <- PCAtest::PCAtest(soil_prop_num, varcorr = T, plot = F, counter = F)
Sampling bootstrap replicates... Please wait
Calculating confidence intervals of empirical statistics... Please wait
Sampling random permutations... Please wait
Comparing empirical statistics with their null distributions... Please wait
========================================================
Test of PCA significance: 17 variables, 44 observations
1000 bootstrap replicates, 1000 random permutations
========================================================
Empirical Psi = 55.8085, Max null Psi = 9.1059, Min null Psi = 4.3334, p-value = 0
Empirical Phi = 0.4530, Max null Phi = 0.1830, Min null Phi = 0.1262, p-value = 0
Empirical eigenvalue #1 = 7.28487, Max null eigenvalue = 3.01024, p-value = 0
Empirical eigenvalue #2 = 3.15651, Max null eigenvalue = 2.34796, p-value = 0
Empirical eigenvalue #3 = 2.59627, Max null eigenvalue = 2.04867, p-value = 0
Empirical eigenvalue #4 = 1.16046, Max null eigenvalue = 1.85095, p-value = 1
Empirical eigenvalue #5 = 0.93253, Max null eigenvalue = 1.63065, p-value = 1
Empirical eigenvalue #6 = 0.74202, Max null eigenvalue = 1.51873, p-value = 1
Empirical eigenvalue #7 = 0.39204, Max null eigenvalue = 1.30991, p-value = 1
Empirical eigenvalue #8 = 0.24218, Max null eigenvalue = 1.18896, p-value = 1
Empirical eigenvalue #9 = 0.19075, Max null eigenvalue = 1.04599, p-value = 1
Empirical eigenvalue #10 = 0.08962, Max null eigenvalue = 0.95824, p-value = 1
Empirical eigenvalue #11 = 0.07595, Max null eigenvalue = 0.83144, p-value = 1
Empirical eigenvalue #12 = 0.05693, Max null eigenvalue = 0.76089, p-value = 1
Empirical eigenvalue #13 = 0.0369, Max null eigenvalue = 0.65485, p-value = 1
Empirical eigenvalue #14 = 0.02507, Max null eigenvalue = 0.56529, p-value = 1
Empirical eigenvalue #15 = 0.01767, Max null eigenvalue = 0.51694, p-value = 1
Empirical eigenvalue #16 = 0.00018, Max null eigenvalue = 0.3948, p-value = 1
Empirical eigenvalue #17 = 5e-05, Max null eigenvalue = 0.30934, p-value = 1
PC 1 is significant and accounts for 42.9% (95%-CI:39.1-49.6) of the total variation
PC 2 is significant and accounts for 18.6% (95%-CI:16.7-24.4) of the total variation
PC 3 is significant and accounts for 15.3% (95%-CI:10.2-17.9) of the total variation
The first 3 PC axes are significant and account for 76.7% of the total variation
Variables 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, and 17 have significant loadings on PC 1
Variables 4, 8, 10, 12, and 14 have significant loadings on PC 2
Variables 5, 6, and 7 have significant loadings on PC 3
Variables 8, 10, 13, 15, and 17 have significant correlations with PC 1
Variables 8, 10, 12, and 14 have significant correlations with PC 2
Variables , and have significant correlations with PC 3Code
pval <- c()
for (i in seq_len(length(pca_test_res$`Percentage of variation of empirical PC's`))) {
obs <- pca_test_res$`Percentage of variation of empirical PC's`[i]
null_model <- pca_test_res$`Percentage of variation of randomized data`[, i]
pval[i] <- (sum(obs < null_model) + 1) / (1 + nrow(pca_test_res$`Percentage of variation of randomized data`))
}
pca_test_pval_adj <- p.adjust(pval, method = "BH")
pca_test_pval_adj [1] 0.005661006 0.005661006 0.005661006 1.000000000 1.000000000 1.000000000
[7] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
[13] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000Code
res_pca <- dudi.pca(soil_prop_num, scannf = F, nf = 4)Code
p_pca_variable <- fviz_pca_var(res_pca,
col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
p_pca_variable_axe3_4 <- fviz_pca_var(res_pca,
col.var = "cos2", axes = c(3, 4),
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
res_pca_var <- get_pca_var(res_pca)
# Fig. 2
(p_pca_variable + p_pca_variable_axe3_4) / free((factoextra::fviz_eig(res_pca)) +
(ggcorrplot::ggcorrplot(res_pca_var$cor)))
Code
practice <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[5]
}))
terroir <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[2]
}))
p_pca_ind_terroir <- fviz_pca_ind(res_pca,
habillage = terroir,
label = FALSE,
repel = TRUE,
mean.point = FALSE,
addEllipses = TRUE,
pointsize = 3,
ellipse.type = "convex"
) +
scale_shape_manual(values = rep(16, length(unique(terroir))))
p_pca_ind_terroir$data$terroir <- terroir
p_pca_ind_terroir$data$practice <- practice
p_pca_ind_terroir_practice <- p_pca_ind_terroir + geom_point(aes(x = x, y = y),
shape = case_when(
practice == "Organic" ~ 2,
practice == "Conventional" ~ 6,
practice == "Conversion" ~ 5
), size = 3.5,
color = rgb(0, 0, 0, 0.7)
)
# value
kw1_terroir <- kruskal.test(res_pca$tab[, 1], terroir)
kw2_terroir <- kruskal.test(res_pca$tab[, 2], terroir)
kw3_terroir <- kruskal.test(res_pca$tab[, 3], terroir)
kw1_practice <- kruskal.test(res_pca$tab[, 1], practice)
kw2_practice <- kruskal.test(res_pca$tab[, 2], practice)
kw3_practice <- kruskal.test(res_pca$tab[, 3], practice)
# Figure 3
p_pca_ind_terroir_practice +
annotate("text", x = 2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - terroir): pval=", round(as.numeric(kw1_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - terroir): pval=", round(as.numeric(kw2_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4, label = paste("Kruskal-Wallis (axe3 - terroir): pval=", round(as.numeric(kw3_terroir$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - practice): pval=", round(as.numeric(kw1_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - practice): pval=", round(as.numeric(kw2_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4, label = paste("Kruskal-Wallis (axe3 - practice): pval=", round(as.numeric(kw3_practice$p.value),4)), size = 3)Warning: Removed 21 rows containing missing values or values outside the scale range
(`geom_point()`).
Code
# Figure S1
ggplot(d_rs@sam_data, aes(x = practice, y = as.numeric(Cu), fill = practice)) +
geom_violin() +
geom_jitter()Warning: Removed 23 rows containing non-finite outside the scale range
(`stat_ydensity()`).Warning: Removed 23 rows containing missing values or values outside the scale range
(`geom_point()`).
Code
kruskal.test(as.numeric(d_rs@sam_data$Cu), d_rs@sam_data$practice)
Kruskal-Wallis rank sum test
data: as.numeric(d_rs@sam_data$Cu) and d_rs@sam_data$practice
Kruskal-Wallis chi-squared = 0.72149, df = 2, p-value = 0.6972Code
d_with_pca <-
clean_pq(subset_samples_pq(d_rs, d_rs@sam_data$paired_name %in% rownames(soil_prop_num)))
res_pca_ind <- get_pca_ind(res_pca)
if (!identical(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name), sort(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name)))) {
stop("ERROR")
}
d_with_pca <- add_info_to_sam_data(d_with_pca, res_pca_ind$coord)Code
dim_df <- as_tibble(d_with_pca@sam_data) |>
select(c(starts_with("Dim"), "nb_seq", "nb_otu"))
cor_results <- correlation::correlation(dim_df)
cor_results %>%
summary(redundant = TRUE) %>%
plot()
Code
lon_lat_with_pca <- d_with_pca@sam_data[, c("long", "lat")]
lon_lat_with_pca$long <- as.numeric(lon_lat_with_pca$long)
lon_lat_with_pca$lat <- as.numeric(lon_lat_with_pca$lat)
MEM_with_pca <- dbmem(lon_lat_with_pca, MEM.autocor = "non-null")
d_with_pca@sam_data$MEM_1 <- MEM_with_pca[, 1]
d_with_pca@sam_data$MEM_2 <- MEM_with_pca[, 2]Code
res_ado_spatial_soil_bray <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE
)
res_ado_spatial_soil_robAit <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_soil_bray_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE
)
res_ado_spatial_soil_robAit_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.01739 0.0086972 0.8898 0.4152
Residuals 72 0.70378 0.0097748 Code
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.23494 0.018072 1.4539 0.1618
Residuals 61 0.75823 0.012430 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 66.0 32.985 0.4173 0.6604
Residuals 72 5691.6 79.050 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 2471.66 190.127 13.793 1.288e-13 ***
Residuals 61 840.86 13.785
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# TABLE 2
kbl(res_ado_spatial_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 0.7692961 | 0.0353392 | 3.462441 | 0.001 |
| MEM_1 | 1 | 1.5042988 | 0.0691030 | 6.770535 | 0.001 |
| MEM_2 | 1 | 1.1007448 | 0.0505650 | 4.954223 | 0.001 |
| practice | 2 | 0.8165301 | 0.0375090 | 1.837516 | 0.007 |
| inter_rank | 3 | 0.7332177 | 0.0336819 | 1.100020 | 0.285 |
| rank | 2 | 0.6416712 | 0.0294765 | 1.444014 | 0.060 |
| terroir | 13 | 4.8718242 | 0.2237972 | 1.686697 | 0.001 |
| Residual | 51 | 11.3313406 | 0.5205283 | NA | NA |
| Total | 74 | 21.7689236 | 1.0000000 | NA | NA |
Code
kbl(res_ado_spatial_soil_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 0.7801743 | 0.0603883 | 3.5267019 | 0.001 |
| MEM_1 | 1 | 0.4296466 | 0.0332562 | 1.9421759 | 0.012 |
| MEM_2 | 1 | 1.4246560 | 0.1102734 | 6.4400191 | 0.001 |
| Dim.1 | 1 | 0.6061032 | 0.0469145 | 2.7398304 | 0.003 |
| Dim.2 | 1 | 0.4590091 | 0.0355289 | 2.0749062 | 0.007 |
| Dim.3 | 1 | 0.2666740 | 0.0206415 | 1.2054739 | 0.215 |
| practice | 2 | 0.5569904 | 0.0431130 | 1.2589105 | 0.139 |
| inter_rank | 2 | 0.3033908 | 0.0234835 | 0.6857243 | 0.914 |
| rank | 2 | 0.5008934 | 0.0387709 | 1.1321201 | 0.285 |
| terroir | 11 | 3.1673809 | 0.2451665 | 1.3016216 | 0.016 |
| Residual | 20 | 4.4243845 | 0.3424631 | NA | NA |
| Total | 43 | 12.9193032 | 1.0000000 | NA | NA |
Code
# TABLE S2
kbl(res_ado_spatial_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 1.5265059 | 0.0715964 | 6.9137648 | 0.001 |
| MEM_2 | 1 | 1.3436557 | 0.0630203 | 6.0856097 | 0.001 |
| practice | 2 | 0.9547461 | 0.0447796 | 2.1620913 | 0.003 |
| inter_rank | 3 | 0.6228395 | 0.0292125 | 0.9403099 | 0.549 |
| rank | 2 | 0.6161566 | 0.0288991 | 1.3953309 | 0.082 |
| terroir | 13 | 4.7758878 | 0.2239993 | 1.6638987 | 0.001 |
| Residual | 52 | 11.4811988 | 0.5384928 | NA | NA |
| Total | 74 | 21.3209903 | 1.0000000 | NA | NA |
Code
# TABLE S3
kbl(res_ado_spatial_soil_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 0.6248397 | 0.0496723 | 2.877821 | 0.001 |
| MEM_2 | 1 | 1.4438602 | 0.1147813 | 6.649979 | 0.001 |
| Dim.1 | 1 | 0.6499372 | 0.0516675 | 2.993412 | 0.003 |
| Dim.2 | 1 | 0.5208441 | 0.0414051 | 2.398849 | 0.005 |
| Dim.3 | 1 | 0.3246755 | 0.0258104 | 1.495356 | 0.088 |
| practice | 2 | 0.5707724 | 0.0453742 | 1.314402 | 0.140 |
| inter_rank | 2 | 0.2656581 | 0.0211188 | 0.611770 | 0.966 |
| rank | 2 | 0.5184269 | 0.0412129 | 1.193858 | 0.194 |
| terroir | 11 | 3.1006425 | 0.2464891 | 1.298238 | 0.030 |
| Residual | 21 | 4.5595725 | 0.3624684 | NA | NA |
| Total | 43 | 12.5792291 | 1.0000000 | NA | NA |
Code
# TABLE S4
kbl(res_ado_spatial_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 305.5208 | 0.0226202 | 1.9770369 | 0.001 |
| MEM_2 | 1 | 277.9747 | 0.0205808 | 1.7987853 | 0.005 |
| practice | 2 | 428.4877 | 0.0317245 | 1.3863802 | 0.002 |
| inter_rank | 3 | 417.9781 | 0.0309464 | 0.9015842 | 0.804 |
| rank | 2 | 433.7451 | 0.0321138 | 1.4033908 | 0.010 |
| terroir | 13 | 3607.0130 | 0.2670571 | 1.7954709 | 0.001 |
| Residual | 52 | 8035.8037 | 0.5949572 | NA | NA |
| Total | 74 | 13506.5231 | 1.0000000 | NA | NA |
Code
# TABLE S5
kbl(res_ado_spatial_soil_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 433.9694 | 0.0317449 | 1.4854760 | 0.001 |
| MEM_2 | 1 | 422.1356 | 0.0308793 | 1.4449687 | 0.002 |
| Dim.1 | 1 | 582.6894 | 0.0426238 | 1.9945441 | 0.001 |
| Dim.2 | 1 | 392.5284 | 0.0287135 | 1.3436235 | 0.018 |
| Dim.3 | 1 | 402.5363 | 0.0294456 | 1.3778805 | 0.021 |
| practice | 2 | 604.4120 | 0.0442128 | 1.0344502 | 0.330 |
| inter_rank | 2 | 444.2780 | 0.0324990 | 0.7603811 | 0.954 |
| rank | 2 | 591.4238 | 0.0432627 | 1.0122210 | 0.424 |
| terroir | 11 | 3661.5622 | 0.2678439 | 1.1394106 | 0.022 |
| Residual | 21 | 6134.9747 | 0.4487744 | NA | NA |
| Total | 43 | 13670.5097 | 1.0000000 | NA | NA |
Variance partitionning
Code
# Figure 8a
res_varpart_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
# Figure 8b
res_varpart_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
res_varpart_robAit_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
# Figure S12b
res_varpart_robAit_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Differential abundance analysis
Indicspecies
Terroir
Code
# Figure 10
res_mpt_terroir <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_terroir_df <- res_mpt_terroir$sign
res_mpt_terroir_df$p.adj <- p.adjust(res_mpt_terroir_df$p.value, method = "BH")
res_mpt_terroir_df$OTU_names <- rownames(res_mpt_terroir_df)
res_mpt_terroir_df_signif <-
res_mpt_terroir_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
tax_tab <- as.data.frame(d_rs@tax_table)
tax_tab$otu <- rownames(tax_tab)
res_mpt_terroir_df_signif_taxo <- left_join(res_mpt_terroir_df_signif,tax_tab, by = join_by("OTU_names" == "otu"))
ggplot(
res_mpt_terroir_df_signif_taxo,
aes(
x = OTU_names,
y = name,
alpha = value,
color = Genus
)
) +
geom_point(size = 6) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))+
scale_alpha(range = c(0, 1))
Code
res_mpt_terroir_rg <-
multipatt(as.matrix(rarefy_even_depth(d_rs, rngseed = 626)@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_terroir_rg_df <- res_mpt_terroir_rg$sign
res_mpt_terroir_rg_df$p.adj <- p.adjust(res_mpt_terroir_rg_df$p.value, method = "BH")
res_mpt_terroir_rg_df$OTU_names <- rownames(res_mpt_terroir_rg_df)
res_mpt_terroir_rg_df_signif <-
res_mpt_terroir_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
ggplot(
res_mpt_terroir_rg_df_signif,
aes(
x = OTU_names,
y = name,
alpha = value,
color = stat
)
) +
geom_point(size = 4) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))
practice
Code
res_mpt_practice <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_practice_df <- res_mpt_practice$sign
res_mpt_practice_df$p.adj <- p.adjust(res_mpt_practice_df$p.value, method = "BH")
res_mpt_practice_df$OTU_names <- rownames(res_mpt_practice_df)
res_mpt_practice_df_signif <-
res_mpt_practice_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>Code
res_mpt_practice_rg <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_practice_rg_df <- res_mpt_practice_rg$sign
res_mpt_practice_rg_df$p.adj <- p.adjust(res_mpt_practice_rg_df$p.value, method = "BH")
res_mpt_practice_rg_df$OTU_names <- rownames(res_mpt_practice_rg_df)
res_mpt_practice_rg_df_signif <-
res_mpt_practice_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_rg_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>Session information
Code
sessionInfo()R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
[5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
[7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Paris
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] geodist_0.0.8 ComplexUpset_1.3.3
[3] microbiomeMarker_1.8.0 formattable_0.2.1
[5] kableExtra_1.4.0 adegraphics_1.0-21
[7] ade4_1.7-22 adespatial_0.3-23
[9] sf_1.0-16 indicspecies_1.7.14
[11] DESeq2_1.42.1 SummarizedExperiment_1.32.0
[13] Biobase_2.62.0 MatrixGenerics_1.14.0
[15] matrixStats_1.3.0 GenomicRanges_1.54.1
[17] GenomeInfoDb_1.38.8 IRanges_2.36.0
[19] S4Vectors_0.40.2 BiocGenerics_0.48.1
[21] factoextra_1.0.7 FactoMineR_2.11
[23] vegan_2.6-6.1 lattice_0.22-6
[25] permute_0.9-7 patchwork_1.2.0
[27] knitr_1.46 DT_0.33
[29] metacoder_0.3.7 targets_1.7.0
[31] MiscMetabar_0.9.2 purrr_1.0.2
[33] dplyr_1.1.4 dada2_1.30.0
[35] Rcpp_1.0.12 ggplot2_3.5.1
[37] phyloseq_1.46.0
loaded via a namespace (and not attached):
[1] igraph_2.0.3 mia_1.10.0
[3] Formula_1.2-5 scater_1.30.1
[5] rematch2_2.1.2 zlibbioc_1.48.2
[7] tidyselect_1.2.1 bit_4.0.5
[9] doParallel_1.0.17 BWStest_0.2.3
[11] clue_0.3-65 rjson_0.2.21
[13] blob_1.2.4 stringr_1.5.1
[15] rngtools_1.5.2 S4Arrays_1.2.1
[17] parallel_4.3.3 RNeXML_2.4.11
[19] correlation_0.8.4 png_0.1-8
[21] cli_3.6.3 ggplotify_0.1.2
[23] CVXR_1.0-12 bayestestR_0.13.2
[25] multtest_2.58.0 MultiAssayExperiment_1.28.0
[27] bluster_1.12.0 BiocNeighbors_1.20.2
[29] TreeSummarizedExperiment_2.10.0 adegenet_2.1.10
[31] mime_0.12 evaluate_0.23
[33] tidytree_0.4.6 ComplexHeatmap_2.18.0
[35] ggh4x_0.2.8 stringi_1.8.4
[37] backports_1.4.1 PMCMRplus_1.9.10
[39] lmerTest_3.1-3 gsl_2.1-8
[41] XML_3.99-0.16.1 Exact_3.2
[43] ggVennDiagram_1.5.2 httpuv_1.6.15
[45] Wrench_1.20.0 paletteer_1.6.0
[47] magrittr_2.0.3 leaps_3.1
[49] splines_4.3.3 jpeg_0.1-10
[51] doRNG_1.8.6 wk_0.9.1
[53] rootSolve_1.8.2.4 ggbeeswarm_0.7.2
[55] DBI_1.2.2 withr_3.0.0
[57] class_7.3-22 systemfonts_1.0.6
[59] htmlwidgets_1.6.4 fs_1.6.4
[61] SuppDists_1.1-9.7 SingleCellExperiment_1.24.0
[63] ggrepel_0.9.5 labeling_0.4.3
[65] SparseArray_1.2.4 cellranger_1.1.0
[67] flashClust_1.01-2 rncl_0.8.7
[69] see_0.8.4 lmom_3.0
[71] effectsize_0.8.8 zoo_1.8-12
[73] XVector_0.42.0 decontam_1.22.0
[75] secretbase_0.5.0 foreach_1.5.2
[77] fansi_1.0.6 caTools_1.18.2
[79] grid_4.3.3 data.table_1.15.4
[81] ggtree_3.10.1 rhdf5_2.46.1
[83] irlba_2.3.5.1 gridGraphics_0.5-1
[85] DescTools_0.99.54 base64url_1.4
[87] lazyeval_0.2.2 yaml_2.3.8
[89] survival_3.6-4 crayon_1.5.3
[91] RColorBrewer_1.1-3 tidyr_1.3.1
[93] later_1.3.2 codetools_0.2-19
[95] base64enc_0.1-3 GlobalOptions_0.1.2
[97] shape_1.4.6.1 limma_3.58.1
[99] estimability_1.5.1 Rsamtools_2.18.0
[101] ggside_0.3.1 foreign_0.8-86
[103] pkgconfig_2.0.3 ggpubr_0.6.0
[105] xml2_1.3.6 GenomicAlignments_1.38.2
[107] aplot_0.2.2 scatterplot3d_0.3-44
[109] ape_5.8 viridisLite_0.4.2
[111] xtable_1.8-4 interp_1.1-6
[113] car_3.1-2 highr_0.10
[115] hwriter_1.3.2.1 plyr_1.8.9
[117] httr_1.4.7 rbibutils_2.2.16
[119] tools_4.3.3 broom_1.0.5
[121] beeswarm_0.4.0 htmlTable_2.4.2
[123] checkmate_2.3.1 nlme_3.1-164
[125] PCAtest_0.0.1 lme4_1.1-35.3
[127] digest_0.6.36 numDeriv_2016.8-1.1
[129] adephylo_1.1-16 Matrix_1.6-5
[131] farver_2.1.2 reshape2_1.4.4
[133] yulab.utils_0.1.4 viridis_0.6.5
[135] DirichletMultinomial_1.44.0 rpart_4.1.23
[137] glue_1.7.0 cachem_1.0.8
[139] Hmisc_5.1-2 generics_0.1.3
[141] Biostrings_2.70.3 classInt_0.4-10
[143] mvtnorm_1.2-4 spdep_1.3-3
[145] metagenomeSeq_1.43.0 statmod_1.5.0
[147] biomformat_1.30.0 ScaledMatrix_1.10.0
[149] carData_3.0-5 minqa_1.2.6
[151] ANCOMBC_2.4.0 glmnet_4.1-8
[153] utf8_1.2.4 seqinr_4.2-36
[155] gtools_3.9.5 readxl_1.4.3
[157] datawizard_0.10.0 ggsignif_0.6.4
[159] gridExtra_2.3 shiny_1.8.1.1
[161] GenomeInfoDbData_1.2.11 energy_1.7-11
[163] rhdf5filters_1.14.1 parameters_0.21.6
[165] RCurl_1.98-1.14 memoise_2.0.1
[167] rmarkdown_2.26 scales_1.3.0
[169] gld_2.6.6 svglite_2.1.3
[171] phylobase_0.8.12 rstudioapi_0.16.0
[173] cluster_2.1.6 hms_1.1.3
[175] ggcorrplot_0.1.4.1 munsell_0.5.1
[177] colorspace_2.1-0 rlang_1.1.4
[179] s2_1.1.6 DelayedMatrixStats_1.24.0
[181] sparseMatrixStats_1.14.0 circlize_0.4.16
[183] scuttle_1.12.0 mgcv_1.9-1
[185] xfun_0.43 multcompView_0.1-10
[187] coda_0.19-4.1 e1071_1.7-14
[189] TH.data_1.1-2 plotROC_2.3.1
[191] iterators_1.0.14 emmeans_1.10.1
[193] abind_1.4-5 tibble_3.2.1
[195] treeio_1.26.0 gmp_0.7-4
[197] Rhdf5lib_1.24.2 DECIPHER_2.30.0
[199] bitops_1.0-7 Rdpack_2.6
[201] ps_1.7.6 promises_1.3.0
[203] RSQLite_2.3.6 sandwich_3.1-0
[205] DelayedArray_0.28.0 proxy_0.4-27
[207] Rmpfr_0.9-5 compiler_4.3.3
[209] statsExpressions_1.5.4 prettyunits_1.2.0
[211] boot_1.3-30 beachmat_2.18.1
[213] BiocSingular_1.18.0 units_0.8-5
[215] MASS_7.3-60.0.1 kSamples_1.2-10
[217] progress_1.2.3 uuid_1.2-0
[219] BiocParallel_1.36.0 insight_0.19.11
[221] R6_2.5.1 rstatix_0.7.2
[223] fastmap_1.1.1 multcomp_1.4-25
[225] vipor_0.4.7 rsvd_1.0.5
[227] ggstatsplot_0.12.3 nnet_7.3-19
[229] gtable_0.3.5 KernSmooth_2.23-22
[231] latticeExtra_0.6-30 deldir_2.0-4
[233] htmltools_0.5.8.1 RcppParallel_5.1.7
[235] bit64_4.0.5 lifecycle_1.0.4
[237] processx_3.8.4 spData_2.3.0
[239] nloptr_2.0.3 callr_3.7.6
[241] vctrs_0.6.5 zeallot_0.1.0
[243] ggfun_0.1.4 sp_2.1-4
[245] pillar_1.9.0 gplots_3.1.3.1
[247] prismatic_1.1.2 ShortRead_1.60.0
[249] locfit_1.5-9.9 jsonlite_1.8.8
[251] expm_0.999-9 GetoptLong_1.0.5